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Abstract 

The dark matter in the universe can be in the form of a superheavy matter species (wimpzilla). 
Several mechanisms have been proposed for the production of wimpzilla particles during or im- 
mediately following the inflationary epoch. Perhaps the most attractive mechanism is through 
gravitational particle production, where particles are produced simply as a result of the expansion 
of the universe. In this paper we present a detailed numerical calculation of wimpzilla gravi- 
tational production in hybrid-inflation models and natural-inflation models. Generalizing these 
findings, we also explore the dependence of the gravitational production mechanism on various 
models of inflation. We show that superheavy dark matter production seems to be robust, with 
Qxh"^ ~ (Mx/10"GeV)2(rRi^/109GeV), so long as Mx < Hi, where Mx is the wimpzilla mass, 
Trh is the reheat temperature, and Hj is the expansion rate of the universe during inflation. 



PACS numbers: 98.80.Cq, 95.35.+d, 4.62. +v; FERMILAB-Pub-Ol/047-A; MCTP-01-16; |hep-ph/010410C 
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I. INTRODUCTION 



The case for dark, nonbaryonic matter in the universe is today stronger than ever [|T|. 
The observed large-scale structure suggests that dark matter (DM) accounts for at least 30% 
of the critical mass density of the universe pc = SH^Mpi/Sn = 1.88 x 10^^^ g cm~'^, where 
Hq = lOO/i kmsec"^ Mpc~^ is the present Hubble constant and Mpi is the Planck mass. 

Despite this compelling evidence, the nature of the DM is still unknown. Some funda- 
mental physics beyond the Standard Model (SM) is certainly required to account for the cold 
and slowly moving particles, X, composing the the bulk of the nonbaryonic dark matter. 

The most familiar assumption is that dark matter is a thermal relic, i.e., it was initially 
in chemical equilibrium in the early universe. A particle species, X, tracks its equilibrium 
abundance as long as reactions which keep the species in chemical equilibrium can proceed 
on a timescale more rapid than the expansion rate of the universe, H. When the reaction 
rate becomes smaller than the expansion rate, the particle species can no longer track its 
equilibrium value. When this occurs the particle species is said to be "frozen out." The 
more strongly interacting the particle, the longer it stays in local thermal equilibrium and 
the smaller its eventual freeze-out abundance. Conversely, the more weakly interacting 
the particle, the larger its present abundance. If freeze out occurs when the particles X are 
nonrelativistic, the freeze-out value of the particle number per comoving volume Y is related 
to the mass of the particle and its annihilation cross section (here characterized by ctq) by 
y oc {1/MxMpiao) where Mx is the mass of the particle X. Since the contribution to 
= Px/pc is proportional to Mxnx, which in turn is proportional to MxY, the present 
contribution to Qx from a thermal relic roughly is independent of its mass and depends 
only upon the annihilation cross section. The cross section that results in Qxh"^ ~ 1 is of 
order 10"'^'' cm^, which is of the order the weak scale. Many theories beyond the SM, e.g. 
super symmetric theories, have stable particles with weak-scale annihilation cross sections, 
and provide candidate weakly interacting massive particles (wimps). 

The simple assumption that dark matter is a thermal relic limits the maximum mass of 
the DM. The largest possible annihilation cross section is roughly M^'^. This implies that 
very massive wimps would have such a small annihilation cross section that their present 
abundance would be too large. Thus, one expects a maximum mass for a thermal wimp, 
which turns out to be a few hundred TeV p. 
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One should note that the computation of the final abundance of the thermal relics assumes 
that the largest temperature of the universe was larger than the relic mass Mx- The thermal 
history of the universe before the epoch of nucleosynthesis is unknown, and the maximum 
temperature in the radiation-dominated phase, dubbed the reheating temperature (Tjih), 
might have been smaller than the mass of the wimp. In such a case, the dependence of the 
present abundance on the mass and the annihilation cross section differs from familiar results 
because of the new parameter T^h [H- This drastically changes the cosmologically allowed 
parameter space of supersymmetric models and re-establishes SM neutrinos as possible dark 
matter candidates p. 

While a thermal origin for wimps is the most common assumption, it is not the simplest 
possibility. It has been recently pointed out that DM particles might have never experienced 
local chemical equilibrium during the evolution of the universe, and that their mass may be 
in the range 10^^ to 10^^ GeV, much larger than the mass of thermal wimps ^ ^, 
Since these wimps would be much more massive than thermal wimps, such superheavy DM 
particles have been called wimpzillas [Q. 

Since wimpzillas are extremely massive, the challenge lies in creating very few of them. 
Several wimpzilla scenarios have been developed involving production during different 
stages of the evolution of the universe. 

wimpzillas may be created during bubble collisions if infiation is completed through a 



first-order phase transition |T^, rT|; at the preheating stage after the end of infiation with 
masses easily up to the Grand Unified scale of lO^^GeV [|12| or even up to the Planck scale 



||13| ; or during the reheating stage after infiation [H with masses which may be as large as 
2 X 10'^ times the reheat temperature. 

WIMPZILLAS may also be generated in the transition between an infiationary and a matter- 
dominated (or radiation-dominated) universe due to the "nonadiabatic" expansion of the 
background spacetime acting on the vacuum quantum fiuctuations. This mechanism was 
studied in details in Refs. Q in the case of chaotic infiation. The distinguishing feature 
of this mechanism is the capability of generating particles with mass of the order of the 
infiaton mass (usually much larger than the reheating temperature) even when the particles 
only interact extremely weakly (or not at all) with other particles, and do not couple to the 
infiaton. 

While the results depend weakly on details such as whether the wimpzilla is a fermion 
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or a boson, or whether it is conformally or minimally coupled to gravity, for the most part 
Qx ~ 1 when the mass of the wimpzilla is approximately the order of the inflaton mass. 
Since hybrid inflation models have (at least) two mass scales and more coupling constants 
than chaotic inflation models, it is worthwhile to study WIMPZILLA production in hybrid 
models [ P^ . 

In this paper we study the gravitational production of wimpzillas after the completion 
of a stage of hybrid inflation. The hybrid scenario involves two scalar fields, the inflaton 
field 0, and the symmetry-breaking field a. Models are parameterized by different mass 
scales and couplings for the two fields. During inflation the inflaton field rolls down along 
a flat potential while the field a is stuck at the origin, providing the vacuum energy density 
driving inflation. However, when becomes smaller than a critical value, (pc, both fields roll 
down very quickly towards their present minima, completing the inflationary phase. It is 
exactly during this phase the gravitational generation of wimpzillas may occur. 

If the WIMPZILLAS are produced at the end of inflation, the fraction of the total energy 
density of the universe in wimpzillas today is given by 



where Hj is the expansion rate of the universe at the end of inflation. Here, flnh"^ ~ 
4.31 X 10~^ is the fraction of critical energy density in radiation today, Tq is the present 
temperature of radiation, and nx{te) is the density of X particles at the time when they 
were produced. The present abundance of the nonthermal wimpzillas is, as expected, 
independent of the cross section P, |^, and one can easily verify that if there is some way 
to create WIMPZILLAS in the correct abundance to give Qx ~ 1, nonequilibrium during the 
evolution of the universe is automatic. 

The paper is organized as follows. In Section II we present some details of the simplest 
hybrid inflation model and discuss the allowed range of the various parameters. In Section 
HI we present our analytical results for wimpzilla production, making use of some general 
results presented in the appendix. Section IV contains our numerical results. Finally, in 
Section V we present our conclusions. 




(1) 
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II. THE HYBRID INFLATION MODEL 



For our computation of WIMPZILLA production, we take the simplest hybrid inflation 
potential as suggested by Linde ||15i [27] 



y (0, a) = ^ {ml - \a^) ' + ^-m^ + l/^V^ . (2) 

This potential has a valley of minima at cr = for > 0c = fna/g- Most of inflation occurs 
while is slowly rolling down from its initial value to 0c. 

During inflation a has a minimum at cr = and its kinetic energy is quickly damped by 
the Hubble expansion. Hence, classically in this naive picture, a remains at for a long 
time before it falls due to some infinitesimal residual displacement of a and/or & about 
[28]. However, this picture is valid, strictly speaking, only when one neglects quantum 
fluctuations. Physically, what will occur is that the quantum fluctuations will grow and the 
long wavelength modes will condense such that different regions of spacetime will behave as 
if they had a classical scalar fleld value of a = ±mo-/ -\/A with domain walls between the plus 
and minus regions. (In the case that the scalar fleld is complex, a cosmic string will form 
instead of a domain wall.) This phenomenon is sometime called spinodal decomposition. 

A relevant observation for gravitational particle production is that the effective stress 
caused by the fleld gradients will increase the pressure of the universe such that the Hubble 
expansion will slow faster. One way to see this is to note that the energy conservation 
equation 

d (pa^) = -Pd{a^) (3) 



tells us that 

' Oi 'l'^ 1 [ dl'^ 

P = Pi — - Pd — , 



(4) 



which implies that a positive increase in the pressure will lead to a faster decrease in the 
energy density, causing a faster decrease in H . Of course, even if the universe contains 
inhomogeneities due to these fleld gradients, one can average over the fluctuation to account 
for an effective energy density and pressure. 

One way of accounting for quantum fluctuations has been presented by Ref. [|1^]. There, 
the canonical formalism is used to quantize the fluctuations about a time dependent zero 
mode (T{t): a = a{t) + 5a{x,t). They argue that the long wavelength modes of 5a{x,t) 
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condense such as to form an effectively homogeneous scalar field Sa{t), whose energy con- 
tribution to the stress energy tensor can dominate over the stress energy of the background 
mode a{t) such that the expansion rate a/ a is damped more quickly than one would naively 
expect from accounting for only a"(t). This effectively homogeneous scalar field 6a{t) has an 
initial condition that is fixed by {6cx^{x,t)) in the background of (f{t). It is 

{S^W - §, (5) 

where the exact numerical factor depends on the boundary condition of the quantum fluctu- 
ations (which cannot be zero due to canonical commutation relations), and Hj is the Hubble 
expansion rate during inflation. 

We will implement this result and simulate the condensation 6a and its fall by letting 
a have a nonzero initial condition at the end of inflation with a value of order Hj/2tt and 
letting it fall, instead of having the condensation component fall. To achieve this, we add a 
perturbation potential 

Vp{<P, a) = BH^ (^a - ^ j exp [-C {<f> - <pf] . (6) 

Then, by adjusting B and C we can simulate the condensate da by making cr(t) roll to the 
new minimum instead. We shall, however, not take into account the potential for 6a {t) as 
is done in Ref. [|17|. In detail, if the potential for a is as given in Eq. (||), the potential in 



which 6a falls would be 

where ex = in our case. Comparing this expression with the tree-level effective potential, 
one finds that the potential for a with a slight displacement from a = achieves the 
same dynamics as 6a{t) if A is replaced with 3A. Hence, if we only consider the case 
where a = forever without the quantum fluctuations, our simulated treatment of spinodal 
decomposition will coincide with that of Ref. [|17| with just the reinterpretation of A — > 3A. 



On the other hand, in reality, since a will never precisely be at zero forever even in the 
nonrealistic absence of quantum fluctuations, a better simulated treatment of the spinodal 
decomposition requires further modifications of the potential along the lines of Eq. (|^) with 
a ^ 0. Since we are primarily concerned with order of magnitude accuracy, and since this 
approximation neglects classical wave scattering effects taken into account in Ref. |18|, we 
will not account for this effective change in the potential for a. 



Let us be more precise about the order of magnitude of B and C. To displace effectively 
a by Hi/2'K at the end of inflation, we must have 



7^2 



B 



Wg 
A 



M, 



PI 



1 



(8) 



ln[l + (?/VCm,]2 

where we have used the COBE determination of curvature perturbations, giving rise to the 
relationship [29] 

A3/2 3.5 X 10-5^3/ ^''> 
Note that the precise value of B and C will not be important to the determination of the 
Bogoliubov coefficient as long as the perturbation potential causes a to fall. We have checked 
this numerically as shown below in the case where we have set m^j = 10~^Mp;. 

We would like to emphasize that while our treatment of spinodal decomposition is ade- 
quate for the purposes at hand, it is far from complete. Since Ref. |T8| argues that generically 
hybrid inflation ends after one oscillation, we cannot realistically probe the parameter space 
in our model where more than one oscillation of the scalar fields is important if we neglect the 
important pressure-related effects due to condensation and classical-wave scattering. Even 
for the one oscillation approximation, the effect of neglecting the pressure due to conden- 
sation and classical-wave scattering underestimates particle production due to the fact that 
the pressure effects increase the nonadiabticity of the expansion of spacetime. Hence, this 
issue certainly deserves more investigation. We note that other related references include 



Refs. M, EO, 21, 221 and references therein. 



The parameters in the potential in Eq. (|^) are constrained by several considerations. 
Constraints on the amplitude and the tilt of the curvature perturbation spectrum generated 
during inflation impose the following constraints on A and g ||16| : 

5 

^ ""-^ ~•3.5xl0-^ (10) 



A3/2 m^M|; 
and 

^— ^ < 0.25. (11) 

The requirement that the cosmological constant term dominates during the inflationary 
regime above 0c imposes a third constraint, 

m\ « (12) 
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Note that the tilt of the curvature perturbation spectrum yields a constraint similar to the 
condition that the field evolution is slow roll; i.e., 

Also, note that the condition that the cosmological constant term dominates during the 
inflationary regime with > 0c also implies the "waterfall" condition (the condition that 
the scalar flelds after reaches 0c roll to the new minima quickly compared to the expansion 
rate). 

With rricr fixed, these constraints collectively determine a region of {g, A) parameter space, 
outside of which is forbidden by the perturbation amplitude and tilt considerations. Yet 
there is one other constraint that we have not discussed. As we have reviewed previously, 
our model does not describe the evolution of the expansion rate of the universe accurately 
beyond one oscillation of the scalar fields after the end of infiation. As we will see in the 
next section, our relic density will depend upon an accurate modeling of the background 
equation for at least one Hubble time at the end of infiation. Hence, our model is valid only 
in the regime in which no more than one oscillation takes place during one Hubble time. 
Let us see how this constrains our parameter space. 

The time scale for the scalar field oscillation is set by the mass matrix (in the (u, 0) basis) 

^ -ml + g'^cj)^ + 3A(j2 2g^(f)a 
^ 2g'^(t)a m^ + ^V, 

which for two extreme values oi a, a = m„/\/\ and cr = 0, becomes 



(14) 



2 r- I 1 2ml ^ 

m (a = mJVA,0 = O) = - 

^ y m; + g^ml/X J 



(15) 



and 

m (7 = 0,0 = o , ■ 16 

We see that although the main oscillation frequency scale is m^, since can be as large as 
0,, = m„/g and since typically m^ -C m^, the actual frequency scale for the oscillation will 
be a weighted time average,whose value can be significantly lower than m^- Let us call this 
weighted average frequency scale m^f, where / < 1 is some constant (typically / is of order 
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10 ^). As far as the Hubble expansion rate at the end of inflation is concerned, in the model 
we study it is given by 



Hi 

Then, the ratio 



\ 3A M|,; ■ VlO-3Mp,y V A 



I-SxIO^Mttt-T^ \^\ GeV. (17) 



fm, V 3A fMpi ^^^^ 
implies that unless rrifj is within a factor of f\f\ of Mp;, many oscillations will occur during 
the one Hubble time when particle production occurs. Hence, the constraint on our param- 
eter space due to limitations of our background field model is that be as close to as Mpi 
as possible. Since Planckian energy densities invalidate semi-classical gravitational physics, 
we will set m„ at the GUT scale, 

= 10-''Mpi, (19) 

assuming that there is some physics separating the GUT scale and the quantum gravity 
scale. Hence, the following interesting set of parameters {g = 0.01, A = 1, = 

^ 652 GeV) which satisfy all the constraints and give a mass scales in the intermediate 
scale (10^^ GeV) and the electroweak scale, cannot be analyzed in our model because in 
this case, Hj/{fmcr) is too small. In fact, even for the single oscillation case, there may be 
some other damping factor for and & which affects the magnitude of H, which of course is 
crucial for the particle production calculation (as we will explain further in the next section). 
Hence, we consider even the numerical calculation results in this article to be only order of 
magnitude accurate. 

Before we map out the parameter space for which our calculation explicitly is valid, we 
would like to show that having nicr close to Mpi forces the scalar fields to have Planck scale 
vevs. This is noteworthy. Because of the possible sensitivity to unknown Planck-suppressed 
operators, scenarios in which the inflaton attains a Planck scale vev may be unattractive 
2§. 



We can model the dynamics of before reaching 0^ as the evolution of a non-interacting 
inflaton in a de Sitter background: 

+ 3if0 + mj0 = 0. (20) 
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For the inflaton field to be slow rolling (overdamped) to the critical value 0c from some 
initial value of (j){t = 0) > (pc we must have m^/H ^1. In that case, taking the least 
damped solution, we have 




(21) 



Note that since 0c = rrifj/g, having ma close to the Planck scale means that 0c will be 
close to the Planck scale. We can be more quantitative by seeing what the constraint 
0(t = 0) < cMpi with c of order unity implies. Since a(t)/a{t = 0) = exp{Hjt), to have 60 
e-folds, we must have 0(t = 0) > 0cexp [60 {m^/ /?>)]. This implies 



where we have taken 0c = nia/g. There are instances when this constraint becomes inde- 
pendent of other constraints. For example, g = 10^^, A = 1, rria = 10"'^ Mp;, = 1.7 x 10^ 
GeV satisfies all other conditions but this one with c = 1. We will neglect this "small field" 
constraint since this is not as fundamental as other constraints. 
In summary, the parameter space that we will explore will be 



The parameter space is shown explicitly in Fig. |l]. 

III. ANALYTIC ESTIMATE OF PARTICLE PRODUCTION 

In the appendix we present a general method of estimating particle production from 
strong gravitational fields. In this section we apply the results from the appendix to the 
hybrid inflationary scenario. 

We show in the appendix that an estimate of particle production requires an estimate of 
the background equation solutions. To start off, let us examine the time variation of H/H. 
After inflation as the scalar fields oscillate about their minima, H/H oscillates. For the 
envelope of the function describing the oscillations we have 




(22) 



3 X 10"^ 



< ^ < 3 X lO^^v^ 



(23) 




(24) 
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FIG. 1: The X-g parameter space in hybrid inflation. The shaded region corresponds to values 



of the parameters allowed by Eqs. (10 - |12|). The lower limit on this region is dashed because it 
represents the "^>" limit in Eq. (p^). 



and the Friedmann equation, 



in 



From these, we find a following general relationship after the end of inflation 



H 



(25) 



(26) 



envelope 

In fact, after the first oscillation the scalar fields will undergo damped oscillation about their 
new minimum, and the scale factor during that time varies in general as 



a{t) = ae(^ 



(27) 



where in the hybrid inflationary case, a ~ 2/3 (which is a typical result of massive scalar 
field oscillation). In reality, this a will have corrections coming from the phase transition 
physics. 

Before inflation ends, the scale factor will be taken to evolve as 



a{t) = aeexp [Hj (t - te)] 



(28) 



with Hi given by Eq. (|T^). 
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Let us now follow the procedure outlined in the appendix to calculate nxite). First, 
consider the contribution to modes that are nonrelativistic at the end of inflation, Ia{k) + 
Ih{k), given in Eq. ( [A8| ) . Assuming Hj is a constant and a{t) evolves as Eq. (^),we find 

(29) 



a. 



4 fc2 e 

where, from the Appendix ti and are defined by 

^physical ('^1 
^physical ('^2. 
Hence, we obtain for Ia{k) the result 



tte a{ti) 
k ap 



-- 2Hi 
Mx. 



Uk) = \ 



1 



M 



X 



2H 



(30) 



(31) 



Next, we calculate the nonrelativistic contribution in the period after inflation, Ib{k), also 
defined in Eq. ([A^). Since ^4. and t^-^-^ are close together, (ts.^^ —t^- <^1/H) for Mx/ H <^ 1, 
and since we are concerned with order of magnitudes, we can just integrate from to 
instead of summing over each i. Since the nonadiabatic region begins at around tg, we take 
^3i ~ te- The final integration time, ^4^, is defined by the condition \H / H\ > Mx- In the 
period after inflation we will take a{t) oc as in Eq. (pTf), so H = a/t and H/H = 1/t. 
Hence, we have 



hik) 



t4» 



a 



H(t) dt = -\n 



«lnf^ 

2 \Mxa 



(32) 



The calculation of the production of modes relativistic at the end of inflation, Ja{k) + Jb{k) 
given in Eq. (|A9|) , is a bit trickier. First of all, consider the contribution Ja{k): 



Ja{k) = \j 



t5 



(33) 



where is the time during inflation when ^physical (^s) ~ '^Hi. During inflation a{t) 
Oecxp [Hi{t - te)], and /i;physical(^5) = ^if/ gives 

k 



ka^exp [-Hi{t^ - t^)] = 2Hj. 



The time is the smallest of the times after inflation when ^physical (^e) - 
^physical (^e) = Mx- In the period after inflation, a(t) = {t/te)^ and H = a/t, so 



MIN 



q/ae 
2Hj/Mx 



l/{a-l) 



1/a 



(34) 
2H or 

(35) 
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where the first term is ^physical '^^s) = and the second term is ^physical ^^^6) = Mx- 
Since will occur after inflation, Ja{k) divides into the parts before and after inflation: 



Ja{k) 



1 

4 

+ - 



dt a^{t) + -\ dt a^{t)H(t) 
2 Jte 



1 - 



M 



X 



2Hi 



2H, 



M 



X 



Of.. 



{2Hj/M^ 



X) 



v2Q/(l-a) 



(g/ae)2/(i-") 



M 



X 



2Hi 



9 



{2Hi 



X 



9 



2Hi 



X 



Or 



(36) 



where ^ is a step function. The second theta function in the second term ensures that 
t^/te > 1- Note that Ja matches la in the limit g/og 1. 

To calculate Jb, we follow the similar procedure as we did for lb, and integrate from 
hi = ^6 to tsiv = ^4jv Note that this is nonzero only when tj-^ < ts^. Hence, we have 



^Inf^VllnfA 

2 \MxaJ 2 \ae 



9 



( Hi 



X 



q_ 



9 



2H, 



M 



X 



(37) 



l/a 



The first 9 comes from tj^ < tg^ and the second 9 comes from using tr^ = tg = {(l/dej 
(see Eq. ^). Note that t = te (o'/'^e)^^" is the time at which the momentum becomes 
nonrelativistic, and it is precisely this regime during which Jb is calculated. If q/ae > 
{2Hj/Mx)°' , then the momentum becomes relativistic and there is no extra contribution 
to Jb- From now on we will assume that a > 1/2, in which case the second 9 function is 
irrelevant. 
Writing 

= (la + h) 0(1 - q/ae + e) + (4 + Jb) e^q/a, - 1) (38) 

where the e indicates that we take the first term when q/a^, = 1, we can finally obtain the 
number density nx{te) through 

2 



nx{te 



Ml 
27r2 

27r2 



1 

48 



1 - 



e 

-1 +—\n'{az) + Ai + A2 + Bi 



(39) 



where z = Mx/Hj and a = 2/3 in our case. The first two terms are the nonrelativistic 
contribution, and the relativistic contributions Ai and Bi are 

(^2 _ 4)2 



A, 



768 



2x 3a 

z. 
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0.08 



0.06 - 0.02z^ + O.OOlz^ 



+ 



(for a = 2/3), 
-24(1 + 3q;)(1 - a)z'^ + (1 + 3q;)(1 - 3a)z^ 



(40) 



\ Iraz 

768(9a2 - 1) 

-48(l-3a)(l-a)] 0^ | 

^ - 0.04 + 0.03;^ - O.OOS;^^ (for a = 2/3), 



(41) 



and 



Si 



1 

54 
0.04 



azj 



3a 



- 1 



+ 



Q!ln(2;) ln(Q;) ln(2;) a'^ln'^{z) Q;ln(Q;) Q;^ln^(a) 



18 6 12 

- 0.04 + 0.07 ln(2) - 0.04 ln2(^) (for a = 2/3). 



+ 



18 



12 



(42) 



We have neglected cross terms as well since we have neglected any phase information 
{i.e., if Pq — Ja + Jb, then was taken to be + J^, which should give a lower bound 
and the correct order of magnitude since both and are positive) . The important result 
is that for small z, one can approximate 

(43) 



In the z < 1 limit, the largest contribution comes from the Ja{k) and Jb{k) terms. This 
corresponds to production of modes that are relativistic at the end of inflation, with approx- 
imately equal contributions to the final value of coming just before and just after the 
end of inflation. We see how the exact behavior of H/H after inflation is important. 
Finally, putting everything together, in the limit z = Mx/Hj <^ 1: 



Trh\ Snznxite 



To y 3 Mj.iHi 



10^) (t^) 

4 ( Trh \ 



2 X 10' 



vio^; VTF^; A (^^^"^ 



(44) 



where the expression is vahd only if Mx < Hi = 1.8 x W^{m^ / IQ-^ Mpif^m-"^ / Xf/^GeV . 
As shown in the next section, the numerical results corroborate this analytic estimate. 

Note that with Mx ~ lO^^Qgy and Trh ~ lO^GeV, we have Qxh^ of order 10"^ 
Characteristic of gravitational production, it is possible to produce dark matter many orders 
of magnitude in excess of Trh- 
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We have left the a dependence in most of the expressions to indicate that the mass 
scahng is sensitive to the fact that the scalar fields enter a regime just after inflation in 
which the scale factor evolves as a matter-dominated universe. The physics of the spinodal 
decomposition is expected to change this effective a, but one would generically expect a 
somewhere between 1/2 and 2/3, which means that the number density of particles pro- 
duced will roughly remain the same. Hence, even though all of our calculations have some 
sensitivity to more than one oscillations (as can be seen in our estimation procedure), as 
long as the scale factor enters a scaling regime at the end of inflation, our results will give 
the correct order of magnitude. 



IV. NUMERICAL CALCULATION OF PARTICLE PRODUCTION 

In this section we describe the results of our numerical analysis of gravitational particle 
production in the hybrid inflation model. The basic hybrid potential was given in Eq. (^. 
As discussed above, the end of inflation is triggered by some perturbation, which we model 
by adding to the basic potential a "perturbed" potential given in Eq. (^). The first issue is 
whether our results are sensitive either to the nature of the end of inflation or the way we 
model it. 

A straightforward exercise is to investigate the sensitivity of particle production to the 
parameters B and C in the perturbed potential we use to trigger the end of inflation. In 
Fig. ^ we show the time evolution of the Bogoliubov coefficient for different choices of B 
and C . As shown in the figure, our results are insensitive to B and C as long as they are 
chosen so as to make Vp negligible outside a very small region around = 0c- Note that we 
also set nifj = 10^^ Mpi for all the numerical work. 

The fact that the final results are insensitive to the exact values of B and C suggests 
(but of course does not guarantee) that gravitational particle production in hybrid inflation 
will be independent of the mechanism that triggers the end of inflation. 

The evolution of the background fields a and determine the expansion rate and the 
change in the expansion rate. Figure ^ is an example of the evolution of the two fields in 
hybrid inflation. For the parameters of this model {g = 0.01, A = 1), the critical value of 
is 0c = O.lMpi. An instability in the trigger field a (driven by the "perturbed" potential) 
causes a and to evolve rapidly to their minima (0 = 0, |cr| = nia-/ V\ = 1) once < 0c. 
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10 20 30 40 

H,(t-t,„,) 

FIG. 2: The absolute square of the Bogohubov coefficient as a function of time for several 
different values of B and C {B is dimensionless and C is in units of Mp^.) We have g = 0.01, 
A = 1, Mx = 0.1-ffo) and k = O.laiHo. The lines correspond to the deformation parameters 
{B,C) given by (1,109); (103,10^); (lO^.lO^); (10^105); (10^,10^); (10^109). The lines are hard 
to distinguish on this scale and asymptotically approach within 10% of each other. 




FIG. 3: An example of the evolution of the inflaton field (p (solid) and a (dashed) as a function 
of time at the end of hybrid inflation. The parameters chosen were ^ = 0.01 and A = 1. 
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Mx/H, 



FIG. 4: The gravitational production of particles during hybrid inflation as a function of A, with 
g set to 0.001. The curves correspond to A as follows: solid, A = 0.001; dots, A = 0.01; dashes, 
A = 0.1; dash-dot, A = 1. The magnitude of O,xh^{TjiH/10^GeV)~^ scales roughly as A~^. 




Mx/H, 

FIG. 5: The gravitational production of particles during hybrid inflation as a function of g, with 
A set to 1. The curves correspond to as follows: solid, g = I (note that this is outside the allowed 
region of g, A parameter space); dots, g = 0.01; and dashes, g = 0.001. 
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Mx/H, 



FIG. 6: gravitational production of particles during natural inflation, with A = 10 ^Mp; and 
U = O.QMpi. 

To calculate the relic density of stable particles produced gravitationally, we integrated 
the background and X-particle mode equations for several different points within the allowed 
regions of parameters shown in Fig. |1], as well as for A = (? = 1, which is well outside it. Our 
results are summarized in Figs. ^ and |^. 

All the curves look similar in form to the mass spectrum for chaotic inflation with a 
potential V{(f)) ~ "^^^^^ The value of Qxh"^ increases with z = Mx/Hj for z < 1, then 
decreases exponentially for z > 1. The reason for this behavior is discussed in this paper for 



the small- 2; region, and in in the large-z limit. 



The numerical results are in qualitative agreement with the result of Eq. (^J)- 
As another example of a single-field model, in Fig. ^ we show the mass spectrum for 
natural inflation ||2^. In natural inflation the potential is usually chosen to be 

V{(l)) = A^ [1 - cos (0/ V2/^) ] . (45) 

Normalizing the parameters to produce the observed temperature fluctuations, a reasonable 
choice of parameters is A = 10~^Mpi and = O.QMpi. With these choices, Hj = 5.1 x 
10-^Mp;. 

As in the hybrid inflation case, in the low- 2; limit VLxh? oc M|^. Again, the numerical 
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results are reasonably represented by Qxh'^ ~ (Mx/10^^GeV)^(T/j///10^GeV). 

V. CONCLUSIONS 

The expansion rate of the universe during inflation, Hj, may signal a new mass scale in 
physics. The particle spectrum of this new mass scale is completely unknown. There may 
be no particles with this new mass scale; an example of such a model is 0^ chaotic inflation. 
There may be only one particle with this mass scale; for example, the inflaton mass in 0^ 
chaotic inflation. Nevertheless, it is very reasonable that one might expect a rich spectrum 
of particles of this mass scale. If this is the case, there may be nearly stable particles of this 
mass scale p6[. Independently of the coupling of the stable particle, they will be produced 



as a result of the expansion of the universe acting on vacuum quantum fluctuations. It was 
shown in Refs. |^, [1^ that such particles would be excellent candidates for dark matter. 
Since the dark-matter particle would have a much larger mass than usual thermal WiMPS, 
they have been named wimpzillas. 

The WIMPZILLA scenario for dark matter seems to be quite robust. The wimpzilla may 
be minimally coupled or conformally coupled, it may be a boson or a fermion, it may couple 
to the inflaton or may be uncoupled. 

The sensitivity of wimpzilla production to the inflation model is one of the subjects 
of this paper. Previous calculations have employed a chaotic inflation model. Here, we 
extend our studies on wimpzilla production to hybrid models and natural-inflation models. 
We have also developed analytic techniques that should provide reasonable estimates for 
wimpzilla production in the limit that Mx < Hj. 

The general picture for wimpzilla production now emerges, and it seems to be rela- 
tively insensitive to the inflation model. The characteristic expansion rate during infla- 
tion, Hj, controls the maximum mass that efficiently can be produced. In all inflation 
models with continuous H, the production of particles with mass larger than Hj is expo- 
nentially suppressed. For particles of mass smaller than Hj, the contribution to Vtxh? is 
(Mx/10iiGeV)2(T«^/109GeV). 

This last expression for Vtxh"^ well illustrates that wimpzilla masses much in excess 
of the reheat temperature may be dark matter. For instance, if T^h = lO^GeV, then 
Mx = lO^^GeV would give Qxh^ in the desirable range. 
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While interesting behavior after inflation like preheating or spinodal decomposition in the 
case of hybrid inflation might change the results, we expect the order of magnitude estimate 
to be correct, and for it to be an underestimate of wimpzilla production. 
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APPENDIX A: ANALYTIC DETERMINATION OF PARTICLE PRODUC- 
TION 



Consider a minimally coupled scalar field with mass Mx ■ The equation of motion for the 
field is 

X + 3HX - ^V^X + M^X = 0, (Al) 

where H is the expansion rate. The scalar field may be expressed in terms of Fourier modes 

= hk/a {a is the scale factor) as 



X 



(27r)3/2a 



ik 



(A2) 



with the usual normalization condition of the creation and annihilation operators, 
S^{k — l), the mode functions hk satisfies the equation 

2 



a^, at 



hk + Hhk 



a \a 



Ml 



(A3) 



In terms of Bogoliubov coefficients and (3k., the mode functions can be written as 

hk = --^==exp (^-i j iOk a~^{t) dt^ + --^== exp (^-i j Uk a'^{t) dtj , (A4) 
where u;l — k'^ + M^a'^. Solving for the mode functions is equivalent to solving the system 



exp ^2i J cuk a ^{t) dtj (5k 



^k 

2uJk 

^k 



(A5) 
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The gravitational production of particles can be expressed in terms of the Bogoliubov coef- 
ficient Pk as 

1 POO A/f"^ 
"-=2A^/o ''"^I'^'l =2^/ ''"^I'^'l • (^*" 
where q = k/Mx = ^physical '^/'^^ with k the comoving momentum and ^physical = ^Z*^ 
the physical momentum. 

The Bogoliubov coefficient to leading adiabatic order can be expressed as (see Ref. ||2^ 
and references therein) 



/■* 1 


\ Hit') 1 






exp {^2iMx J 


1 + gVa2(t'). 



dt"^l + qya^{t")] . (AT) 



This formula breaks down when \j3\ is of order unity (which may occur in our scenario), but 
let us use it to estimate the order of magnitude scales. 

The magnitude of f3q depends on the magnitude of the argument of the exponential in 
Eq. (|A^). If the argument is of order unity or greater, then the oscillatory behavior will 
damp \Pq\. Thus, the final magnitude of Pq depends on the size of q/a{t). This leads to a 
natural division of particle production into the cases where q/a{t) = ^physical /^-^ larger 
or smaller than unity. We will consider the two cases in turn. 

First consider production of nonrelativistic particles: q/a{t) = ^physical/^^ ^ ^' This 
case further splits into two subclasses. 

Mx < \H/H\: In this case, the oscillations are not important, and one simply integrates 
H{t) to the point that it becomes negligible. 

Mx > \H / H\: In this case, the oscillations cancel most of the contribution to the inte- 
grand. 

Now consider production of relativistic particles: q/a{t) = ^physical /^-''^ ^ ^' ■'■^ ^^^^ 
case, the frequency of the oscillations just becomes the physical momentum. Again, this 
case divides into two subclasses. 

^physical ^ I^MI' where A(t) = H{t)/ [1 + q^ /a?{t)]: Since q/a(t) > 1, this is equivalent 
to ^p]2ysical ^ ^H + H/H . In this case the oscillations are not important, and one simply 
integrates H{t)/[q'^ /a^it)] to the point that it becomes negligible. 

^physical ^ + H/H : In this case the oscillations cancel most of the contribution to 
the integral. 

We will neglect the marginal case of q/a = 1, since this will be roughly accounted for in 
the estimates of the above cases. The different cases and subcases are given in Table |. 
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TABLE I: This table summarizes the different cases in the analytic calculation of gravitational 
production of particles. 



relativistic /nonrelativistic 



subcase 



oscillations 



/5 



nonrelativistic 
nonrelativistic 
relativistic 
relativistic 



Mx < \H/H\ 
Mx > \H/H\ 
^physical < + H/H 



none / dtH{t) 

many damped 

none J dtH{t) [q/a{t)y 
many damped 



A key to developing analytic approximations is the behavior of H and \H/H\. During 
inflation, H is roughly constant (denoted as Hj) and H/H ^ Hj. After inflation, H/H is 
negative, and oscillates (with decreasing amplitude) between zero and approximately —H. 
This behavior is illustrated in Fig. |^ in the simple chaotic inflation scenario. During the 
matter-dominated (MD) phase and the radiation-dominated (RD) phase, H/H = —3H and 
H/H = -AH, respectively, so 2H + H/H =[1 (MD) or 2 {KD)]xH 

Since ~ during inflation, from Table | we see that production of nonrelativistic 

particles is suppressed during inflation and production of relativistic particles during inflation 
is suppressed if ^pliysical 

> H. 

Let us now turn to the estimate of the number density. The particular inflation model, 
together with the behavior of the expansion rate immediately after inflation, will determine 
the efficiency of gravitational particle production. Here we will give a recipe that can be 
adapted for several models. 

We are mainly concerned with the case when Mx/Hj < 1. (Particle production is 
exponentially suppressed for Mx/Hj > 1: this case was addressed in detail in Ref. 
The result will depend on whether the particle was relativistic or nonrelativistic at the end 
of inflation. 

First, consider momentum modes where the particle was nonrelativistic at the end of 
inflation, ^physical^^^) — "^^^ calculation divides into production during inflation 

and post-inflation production. During inflation, the growth in is only when Mx < 
^physical ^ 2iJ/. After inflation, the particle is nonrelativistic, and growth occurs during 
periods when Mx < H/H. Using the results summarized in Table | (recall that ^physical ~ 
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FIG. 7: The behavior of H/H and the inflaton field at the end of inflation in a simple chaotic 
inflation model {V ~ M?(/)^). Notice the oscillatory behavior of H/H after inflation. 



k/a = Mxq/a), 

1 rti<t2<U Hit) ^ 1 /■*4i 

A(*physical(* J < M.) - ^ i. * ^ + S 2 * ^(«) 

= 4 + 4. (A8) 

Here, la is the growth during inflation in the interval {^1,^2} where the times are defined 
such that ^physical ('^1) ~ '^^^ ^^"^ ^physical (^2) = Mx- h is the growth after inflation in 
the intervals {tziiti^} when Mx < H/H. 

Note that we have neglected any phase information between various contributions. These 
interfererence terms should be important in only some special cases and not generically 
because in most cases only one term will dominate. 

Now, consider momentum modes where the particle was relativistic at the end of inflation, 
^physical — ^i^^^ mode was relativistic at the end of inflation, it must have 

been relativistic throughout inflation. From Table | we see that during inflation, the growth 
in the amplitude of (3k only occur when 2Hi > ^pjiygjcah After inflation, the mode will 
remain relativistic so long as ^physical ^ ^"^^ continue to grow so long as 

2H > ^physical- After the mode becomes nonrelativistic (^physical ^ Mx) it wih grow 
only during periods when Mx < H/H. Using the results summarized in Table | (recall that 
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^physical = k/a = Mxq/a), 

1 /"^S Hit) 1 f*8 

= + (A9) 

Here, Ja is the growth during and (possibly) after inflation in the interval {tsjte} where the 
times are defined such that ^physical ('^s) ~ ^-^^ smallest of times after infiation 

when either ^physical (^b) = 2if or ^physical (^s) = Mx- h is the growth after infiation in 
the intervals {t^.^t^i^} when the mode is nonrelativistic and Mx < H/H. 

Once again, we have neglected any phase information between various contributions for 
the reason discussed above. 

To use these facts to estimate the relic density produce, one must first obtain a reasonable 
estimate of a{t) from the background equations. 
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